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SUDDEN BENDING OF A CRACKED LAMINATE 
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Institute of Fracture and Solid Mechanics 
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Bethlehem, Pennsylvania 18015 

and 

* 

E. P, Chen 
Sandia Laboratories 
Albuquerque, New Mexico 87115 

ABSTRACT 

♦ 

A number of laminated plate theories have been developed in recent times to 
analyze the static and dynamic response of composite laminates with or without 
the presence of stress concentrators such as holes, cracks, etc. Many of the 
theories tend to quickly become Intractable when considering the determination 
of the state of affairs near the singular crack edges that are present in the 
laminate, particularly if the loading is time dependent. Additional uncertain- 
ties arise due to the lack of information on the mechanical properties of the 
interface through which load transfer takes place between the adjacent layers. 
This paper focuses attention on the Intensification of stresses near a through 
crack In the laminate that suddenly undergoes bending. A dynamic plate theory 
is developed to include many of the essential features of the problem such as 
material nonhomogeneity in the thickness direction, realistic crack edge stress 
singularity and distribution while the parameter dependence of various signifi- 
cant quantities Is also assessed. Of particular interest is the variation of 
the dynamic stress intensity factor with time. Numerical results for different 

*Dr. E. P, Chen was on the faculty at Lehigh University. 
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geometric and material constants are displayed graphically to show how they can 
affect the transfer of load to the vicinity of a through crack In the laminate 
that undergoes sudden bending. 


INTRODUCTION 

The damage of laminated composite materials Is, to say the least, very com- 
plex since It Involves various modes of failure such as fiber breaking, matrix 
cracking. Interface delamination, etc. Analytical modeling would be beyond ap- 
proach if all these failure modes were to be accounted for. The spirit of frac- 
ture mechanics is to assume that a critical single flaw or damage zone exists 
and can lead to Instability in terms of load applied to the laminate, Damage 
accumulated in the composite other than the dominant flaw may often be simulated 
by changing some of the mechanical properties of the composite v^hich are usually 
the stiffness of the constituents. Although not all laminates can be identified 
with a single characteristic damage state, the single-flaw fracture mechanics 
approach will be taken in this analysis in order that a sensitivity study on the 
physical parameters affecting laminate fracture can be made possible. One of 
the main objectives of this investigation is to come forth with a feasible dy- 
namic theory of the laminate plate for analyzing composite failure due to crack 
propagation. 

As a consequence of increased use of laminate composites in aircraft and 
other high speed vehicles, the analysis of the fracture behavior of layered com- 
posites has attracted the attention of a considerable number of investigators 
[1,2]. A variety of diverse approaches has been proposed to analyze laminate 
failure and a collection of papers on this subject can be found in [3], The 
role with which the interfaces play in transferring the load from one layer to 


the next In the laminate was emphasized. Because of the differ'ence in the ma- 
terial properties of the adjacent layers, the stresses across the Interface ex- 
perience steep gradients. Only recently, a comprehensive, study was made on how 
the conditions in the Interface can influence composite failure [4], Even though 
the interface may be relatively thin when compared with other dimensions of the 
composite, the resulting stresses can be sensitive to the material properties 
of the Interface depending on the loading conditions. There exists no theory at 
the present which can relate the strength of a composite structure to the condi- 
tions In the Interface. This aspect of the problem Is emphasized In this re- 
port. 

The aforementioned difficulties become even more overwhelming when the load- 
ing Is time dependent. There Is f<o need to emphasize the virtue for constructing 
approximate dynamic theories for laminate composites, particularly for handling 
crack problems. In the case of bending loads. It Is essential that the three 
physical boundary conditions of bending moment, twisting moment and transverse 
shear stress be satisfied individually on the crack edge. Such a theory has 
been developed by Mindlln [5] for a single layered plate made of isotropic and 
homogeneous materu'"* and applied to solve a number of crack problems [6]. An 
equally effective theory Is described herein for the dynamic bending of laminate 
plates. Each layer of the laminate assumes different elastic properties and 
Is attached to the next layer with continuous strains across the interface. 

The problem of a through crack In a balanced symmetric laminate Is solved for 
a moment applied suddenly on the crack surface. Not only are the qualitative 
features of the three-dimensional stress distribution preserved in the vicinity 
of the crack front, but, perhaps more significantly, the dynamic stress Intensity 
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factor, which is a quantitative measure of the load transmitted to the crack, 
is determined in terms of the significant material and geometric parameters such 
that an effective study an laminate fracture can be made. 


DYN AMIC THEORY OF LAMINATED PLATE 


Without loss in generality, a four layered composite plate will be considered 
as shown in Figure 1, The two middle layers are made of a material with shear 
modulus Poisson's ratio and mass density while the two outer layers 
have the properties p 2 » ^2 ^2* ^ rectangular Cartesian coordinates 

X, y and 2 are attached to the mid-plane of the laminate such that the layer 
•properties are symmetric with respect to the xy-plane with z being the thickness 
coordinate. The total height of the laminate is h with each layer having the 
same thickness h/4. The outer edges of the laminate are sufficiently far away 
from the crack so that their influences can be neglected. 


aw-cc (Uiumptloiu and AeXationi, The layers of the laminate in the thickness 
possess different material properties pj, vj and pj (j * 1,2) such that (p^, 
v^,p^) prevails in the range 0<lz|<h/4 and (p 2 ,V 2 ,P 2 ) applies to h/4<|z|<h/2. 

The surfaces of the laminate are free from tangential tractions 
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'^xz * ’'yz * ^ ^ 


( 1 ) 


but may be subjected to nonnal pressures and qg as follows; 
-qi(x,y,t) for z « h/2 

‘^ z • 

-qg(x,y,t) for z * -h/2 
In the sequel, the notation 


( 2 ) 


q(x.y,t) » q 2 (x,y,t) - q^(x,y,t) 


( 3 ) 


will be used. In plate theory, it is more convenient to work with the moments 
Mx» My, Hj^y and shearing forces Qy defined in the usual manner as 

h/2 


(Mj^,My,Hxy) • ia^fOy,x^y)z(iz 


( 4 ) 


h/2 




-h/2 


xz* yz' 


From the stress and strain relations and equations (4), the expressions 


^ = DiC(r^)^ + v^(ry)^] ^ DgCir^)^ + v2(Fy)^] 
My = Vi(rj^)y] + DgCCPy)^ + 

^xy “ ^ 2 ^ ^l^^xy^^ ^2^^xy^2 


( 5 ) 


and 


Qx = * >'2<'’xzy 

Qy =r 


( 6 ) 
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are developed provided that the quantities (r„) , (r^,) , (r„_) (j ■ 1,2) 

* j > j j 

stand for 


Ofi 

Qc -h/4 

V ^^^y^ * ^ ^‘hjz 

h/2 

* J4 



(7) 


In equations (5), D^ and D 2 are the flexural rigidities of the layers given by 
P,h3 7p,h3 

“1 * 4S(T-v,)’ °2 ■ wrn^ 


The constant k in equation (6) accounts for the thickness-shear motion of the 
plate and takes the value of irZ/TT as given in [5]. 

Now, let the displacements be continuous through the interfaces by letting 

\ = z4>j^(x,y,t), Vy = zij)y(x,y,t), = w(x,y,t) (9) 

Making use of the strain-displacement relations together with equations (7) and 
(9), it is found that 
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3'1'v 9’1'u 

C-.)- * C'x >2 ' ■ V, * 


1 


1 


3'<'v 31'x 


(10) 


(Px^)^ “ ^ '*’ 


8W 

3y 


Hence, the moments My and Hj^y can be expressed in tenns of the displacement 


functions ij»y and w: 


”x ° ‘’o'^siT * '’o 

3V. 3V. 

”y = “oCs/ * '■o 

V * r <’■'’ 0 X 3 / + 

The same applies to and Qy which become 

Ox ' if '•^o(*x ^ I?) 

^y " T? ^’*^Q^’*'y 

Note that Dq, vq and Pq are defined as 

. . . . Di'’1+D2'^2 _ 

Dq = Di + Dg, Vq - , Pq ^ 


( 11 ) 


(12) 


(13) 


Equations (11) and (12) are, in fact, similar to those derived in [5] for the 
case of a single layer homogeneous plate except that the constants D, v and u 
are now replaced by Dq, vq and ug. 

GoveAfiOig cU^^eA.oituit equattona. Consider the elastodynamic equations of 
motion given by 


_JL + JM + 

ax ay az 


a^u. 


a^v. 


llxi + !!i + = p(j) 

ax ay az at^ 


(14) 


3t 3ty 80 - 
TT ^ “3y“ ^ sT • 


a'^Wj 

aF“ 


in which the mass density may vary in the thickness direction of the laminate. 
Multiplying the first two equations by z, expressi?ig the stresses in terms of 
moments and integrating the results with respect to z from -h/2 to h/2 lead to 


DM. 


3H 


ax 

3H 


ay ^x 


ja!. + Q 
ax ay ^y 


- ''0 . . 

- Tf aF~ 


= ^ h3 
12 " aF^ 


(15) 


. 5i 

ax ay 


+ q = ph 


a^w 

IF 


in which 


Pq “ ^ (Pl ^pg) » P = ^ (p*! + Pg) 


. 1 


(16) 
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The result of inserting equations (11) and (12) into equations (15) is a 
system of second order partial differential equations 


^ « 0+Vq) g Bf 3i|> 2* Pq 32,|, 

2-2- + -2-^ Og w * 1?^ “of'l’x * lx> * 17 " aF“ 




"''Q/ „ 3 ^‘^y\ 


7T * 3y' 


PO .afi 

3t^ 


(17) 


“0<’ “ ^ ?r * a/* * “I “ p" Ip- 


where 7 ^ » d^/dx^ + 3^/3y^ is the Laplacian operator in two dimensions. Equa 
tions (17) may be combined to give a single equation 




32 


(“o’^ • "IT'lpX''^ ■ li 


32w 


’^0 


12Dq h2p 2 


(18) 


solving for the transverse displacement w(x,y) of the laminated plate, 

BoundoAtj aonditioifii). In order to derive the boundary conditions that must 
be specified on the crack, consider the energy stored in the laminate 


h/2 

W = / Wdz 

-h/2 


^(M 


x''x^ 


”y'’y * "xy^y ^ «x''xz " Oy’'yz) 


(19) 


in which Ty,-— , 
(10). Equations (5) 


r are related to i|^„ and w as indicated in equations 

A y 

and (6) may thus be applied to render 
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4W » DQ(Uvo)(VFy)2 + |i + rjy] 


( 20 ) 


Since the physical constants Dq(1+Vq) and Dq(I-Vq) are positive, tJ" is a positive 
definite quantity. Hence, W vanishes if and only if the equivalent strains 
Tyt etc., vanish individually. Equation (20) also implies that 


M. 


aw » - iW_ H 
”y ” ar^* "xy 



and 


( 21 ) 


n - n - J}L 
^x ’ ar^^’ ^y ar^^ 


The kinetic energy in the laminate is 



( 22 ) 


which, when expressed in terms of \p^ and w, takes the form 

A y 

T = ^§5- [(— )^ + (■^)^ + (— )^ 


It is now possible to write down the expression for the total energy of the 
laminate at time t; 


_ _ t ^ ^ tTh aw 2 

T + V = / dt // {“ 2 ^ C(glf“) ■*■ (g^) ] 2 - (pf) )dxdy 


+ / dt // 1^ dxdy + Tg + V„ 


(24) 
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where ? Is the total potential energy. Note that Tq and Vg are the values of 
T and ? corresponding to time tg. Equation (24) may be integrated by parts and 
the results may be arranged to read as* 




He 


T*V.i 


‘o ' 


I 

t. 


+ / dt // q II dxdy + Tg + Vg 


(25) 


The above result may be interpreted as the total energy in the laminate at time 
t and consists of the initial energy at tg plus the work done by the external 
forces along the edges and over the surfaces of the laminate during the time in- 
terval t-tg. The initial and boundary conditions for the laminate can now be 
easily extracted from equation (25). They can be summarized as follows: 

(1) On the laminate or crack edges : Any combination containing one member 
of each of the three pairs (^, M^), (-g^, H^g) and (~|, Q^^) may be specified 
on the crack or laminate edge. 


(2) Throughout the laminate: The initial values of ii>^, ij;., and w and their 

A y 

time derivatives need be known. 

(3) Tractions and Displacements : The external load q or the displacement 

w on the laminate may be specified. 

This completes the development of the dynamic laminate plate theory which 
will be used to solve a crack problem. 


*Refer to page 45 for the derivation of equation (25). 
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A CRACKED LAMINATE PLATE 


As an example, consider the laminate In Figure 1 to be initially at rest 
and bent suddenly by a moment with a constant magnitude of Mq maintained on the 
crack surfaces. The conditions can be stated as 

Qy(x,o.t) = Hjjy(x,o,t) » 0 for 0<|x|<« (26) 

and 


My(x,o,t) * -MQH(t) for lx{<a and i|)y(x,o,t) = 0 for Ixl>a (27) 

which Is of the mixed type. The displacement functions are subjected to the 
conditions that 

linj C'1'«(x,y,t), i/»y(x,y,t), w(x,y,t)] = 0 

x2+y2-H» ^ y 

No other external forces or constraints are present. 


Laplace. txan6^oAm. The governing equations (17) will be solved by intro 
ducing the Laplace transform pair 


f*(p) = 7 f(t)e"P^dt 
0 

(28) 

= ST I, 


where the second Integral Is over the Bromwich path. Applying the first of 
equations (28) to (17) yields 
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0-Va) ^ * (Hv«) « 3i|)* 3</»* 

~r^ “olx * 3y^) 

3 /^'^x ®'^y\ 

Tf2h / , * ’^’^x ^'^'y X -.. 2 * 

TT >'o<'' *' *3jr* 9j^> ' '’"P^ 


fr '*ot'''x * l^>' ^ 

^ Pot'i’J f-> ' ^ P^Py 


( 29 ) 


The analysis may be simplified by letting 


★ it 


★ Hr 


(,* . Hu t 9H_, / . »i_ . »tL 

5x 3y * ^y ay ax 


(30) 


such that equations (29) simplify to 


lx 'P'1'* - <P«0 1- s 


■’V* - S"’w*} + -^ly (p2-0j2)H* = 0 


1-Vr 


{72/ - (R«S|^ + s""')/ - S'^’w*} - (V^-(o 2 )H* = 


9y 


(31) 


ir ★ 


V2(<fi +W ) - S'^W = 0 


The new quantities introduced in equations (31) are defined as 


R = hi S = s" = ?• = ^ 

U* ^ ^0 Dq ’ '’*0 Dq 


(32) 


and 


2(R6§fS~V)(D^-fD3) 

(1-Vi )D i+(1-V2)D2 


„2 = 


(33) 
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Furthermore, If 


* (3-1 )w* 

Is Introduced Into equations (31), It can be shown that 
v^w* - o^w* * 0 
while a and 3 are given by 


-1 Sdj; 

a2 • R6g + S + |qp, 3 » ^ 


'fit if 

Consequently, the functions and \py in equations (30) become 


★ 9Wi 3W« 

h * IT* JT ^ 


ii 

By 


* 9W, 3W2 

I'y - (3^-1) jy- + (32"1) ■ 


3 X 


and w* may be written as 

* * ★ 

W = W-| + Wg 

In equations (37), 3i and &2 ^re given as 

8i. 2- a|_, 

In which 

“f,2 ' ? 
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( 34 ) 


(35) 


(36) 


(37) 


(38) 


(39) 


(40) 




( 41 ) 


if ^ it 

It Is now apparent that once H • and W 2 are found from 

(v2-u)2)h ■ 0, (7^-<x|)w^ ■ 0* (72-a|)W2 * 0 

the problem Is basically solved In the Laplace transform plane. 

JntegAal equation. Taking advantage of the symmetry condition with respect 
to the y-axis, It Is not difficult to show that the following Integrals 

* ? * Ml "Yiy 

w,(x»y,p) ■ f / B''^(s,p) cos(sx)e ' ds 

' 0 

wj(x,y,p) - |7 B^^^(s,p) cos(sx)e ^^^ds (42) 

0 

H*(x,y,p) « |-/ B^^^(s,p) s 1 n(sx)e ^^^ds 
^ 0 

satisfy equations (41) provided that 

^1,2 " Y 3 = (s^+(D^)^^^ (43) 

The unknowns B^^^, B^^^ and B^^^ must be determined from the boundary conditions 
in equations (26) and (27) whose Laplace transform are 

Qy(x,o,p) = H*y(x,o,p) = 0 for 0<x<« (44) 

and 

★ ^0 

My(x,o,p) = - p~ for 0<x<a and ijiy(x,o,p) = 0 for x>a (45) 

The appropriate quantities in equations (44) and (45) may be obtained by first 
putting equations (42) into (37) and (38). This gives 
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( 46 ) 


♦x ■ • 1/ (sC(6,-1)B<’>(s,p)e + (e2-1)B(2)(s,p)e'^2^] 

+ Y 3 B^*^s,p)e s1ti(sx)ds 

/ U6,-l)Y,B(’)(s,p)e''''^ + (62-l)YjB(2)(s,p)e'''2^ 

+ sB^^^(s,p)e cos(sx)ds 
and 

W* - f J (B(’)(s.p)e + B<2)(s,p)e'^2^) cos(sx)ds 

The Laplace transform of equations (11) and (12) will clearly Involve (,*, / 
and w*. Equations (46) and (47) and equations (45) can be satisfied If the^ 
function C(s,p) obeys the dual Integral equations 

oo 

/ C(s,p) cos(sx)ds * 0 x>a 
0 “ 

/ sG(s,p)C(s,p) cos(sx)ds = 
with G(s*p) beinq a known function 

0 -vg) 

-J 8(s,p) > ((1 -S,)(y^-V(|S2) /(sy,) - ( 1 -B 2 )(y|-v„sZ)®/(sy 2 ) 
2 sy 3 ( 1 -\>o) (a|-a 2 )/a)^}/(a|-a|) 


( 48 ) 

x<a 
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1 


I 


The conditions In equations (44) may be used to relate the functions 8^^^, 
and to C(s»p).' 


B^^^(s,p) 

B^^^(s,p) 


TjTaf-l^r 

(l-VQ)s2+a| 

- r2{a-fc^) 

sO-voXsj-s,) 

— spsf — 


(50) 


Without going into details, the solution for equations (48) is of the form' [6] 


TfMQa^ 1 


Cis.p) ‘ ^ 


(51) 


where Jq is zero order Bessel function of the first kind and the function T^ls.p) 
can be found from a Fredholm integral equation of the second kind; 

/(5.P) + } L(S,n,p)/(np)dTi * (52) 

0 

The kernel L(5,n.p) is symmetric in ? and n and takes the form 


L(c.n.p) » / sCG(|, p) - l]jQ(s 5 )jQ(sn)ds (53) 

Equation (52) can be evaTu«,ted numerically for 4'*(c,p) in the Laplace transform 
domain and then inverted into the time domain by using the second of equations 
(28). 

Vi/nanUc moment ^temntij ^aeton. The time dependence of the solution may be 


recovered by two different procedures. The first is to apply the Laplace inver 


Sion formula to the quantities of Interest and obtain the complete solution as 
a function of time. Such an approach Is not only cumbersome and can often re> 
suit In a considerable amount of difficulties In numerical calculations. In 
fracture mechanics, since It 1$ only necessary to focus attention on the state 
of affairs near the crack front, Sih et al [7] have suggested to obtain the 
asymptotic stress solution In the Laplace transform domain such that the time 
Inversion Is applied only to the first term of the stress expansion near the 
crack tip. This approach has greatly simplified the analysis and will be used 
here. 


The local solution may be found by expanding the Integral In equation (51) 

fir # 

for C(s,p) for large values of the argument s. Once the moments and 

are expressed In terms of C(s,p), the resulting Integrals may be evaluated to 
give the asymptotic expansions: 


K (d) 

M*(r,0,p) * cos (I - sin |> s1n |^} + 0(r°) 

K^(d) 

M*(r,0,p) * ■■ cos j {1 + sin sin |^} + 0(r°) 

★ 

H*y(r,0,p) » cos I sin | cos |^ + 0(r°) 

fir fir 

where Q and Q are nonsingular and remain finite as r-+0, 1,e., 

X y 


(54) 


Q* = Qy = 0(r°) " (55) 

The polar coordinates r and e are measured from the crack front as shown In Fig- 
ure 1. The parameter 
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Is the Laplace transform of the dynamic moment Intensity factor and f*(l,p) de- 
notes the values of the function f*(5,p) near the crack border c*l. 

Applying the Laplace Inversion theorem to equations (54) yields the solution 
as a function of time: 


K (t) 

Mj^(r,0,t) ■ f f ■*' 

K (t) 

My(r,0,t) * -3 cos I {1 + sin |* sin |^} + 0(r°) 

i^S'r 

K It) 

Hj^y(r»e,t) « — - cos j sin |= cos |^ + O(r^) 


(57) 


The dynamic moment Intensity factor K|(t) may be computed from 



once ^^(l.p) Is known. 

Um&fUcxLt Ke^iittA. Since the procedure for solving the Fredholm Integral 
equation Is already well known, it is not necessary to cover the details. The 

if 

numerical values of S' (l,p) in equation (58) are given In Figures 2 to 4 for 

the three different values of yg/vi = 0-1 > 1-0 ^nd 10.0. The Poisson's ratio 

and mass density for the layers are assumed to be the same as their variations 

In the thickness direction do affect the results appreciably. The function 

f (l,p) is seen to increase monotonically with Cgi/pa where C 21 = (U|/pi) ' 

is the shear wave speed of the material in the outer layers. 
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As an Indication of the load intensity transmitted to the crack edge region 
as a function of time, the moment intensity factor K^(t) will be computed from 
equation (58) by using the results in Figures 2 to 4. Figures 5 tc 7 display 
the variations of the normalized quantity K.|(t)/MQ^ with the dimensionless 
time parameter C 2 -jt/a for M 2 /'^i ' while the crack length to 

laminate thickness ratio 2a/h takes on the values of 1, 2 and 4. Generally 
speaking, K^(t) tends to increase with time reaching a peak and then acquires 
an oscillatory character. The peak value of K-|(t) appears to be inversely pro- 
portional to the ratio of 2a/h, i.e., K^(t) maximum at 2a/h = 1 is larger than 
that at 2a/h * 4. The moment intensity tends to decrease as the crack length 
is increased. Also, K^(t) maximum occurs earlier when the shear moduli in the 
outer layers of the laminate is larger than those in the inner layers. Refer 


to the curves in Figure 7 for ^ ^ those in Figure 5 for U 2 /M 1 < !• 

The influence of can be best illustrated by fixing the ratio of 2a/h and 

use ug/v-i as a varying parameter. Figure 8 shows a plot of K^(t)/MQV? versus 
C 2 -|t/a as takes the values 0.1, 1.0 and 10.0. It is clear that the crack 

edge moment intensity can be reduced by letting pg ^ making the shear 

moduli of the inner layers to be larger than the moduli of the outer layers. 


CONCLUDING REMARKS 

A dynamic laminated plate theory is developed with emphases placed on ob- 
taining effective solution for the crack configuration where the 1//F stress 
singularity and the condition of plane strain are preserved. The radial dis- 
tance r is measured from the crack edge. Although each layer in the laminate is 
assumed to be isotropic, it is a simple extension to include anisotropy simu- 
lating the directional properties of fiber reinforcement. This additional com- 
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plexity was not thought to be necessary in this preliminary analysis. 

Several revealing conclusions can be made from the numerical results of the 
example on the sudden bending of a cracked laminate when compared with a single 
layer homogeneous plate. 

(1) The crack moment intensity tends to decrease as the crack length to 
laminate plate thickness is increased, Hence, a laminated plate has the de> 
sirable feature of stabilizing a through crack as it increases its length at 
constant load. 

(2) The level of the average load intensity transmitted to a through crack 
can be reduced by making the inner layers to be stiffer than the outer layers. 

The foregoing comments are strictly based on the concept of moment intensity 
factor as used in the theory of fracture mechanics. In the normal course of 
design, other considerations must also be accounted for. However, the point has 
been made that the present theory, although approximate, is useful for analyzing 
laminate failure due to crack propagation. 
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CompittcA ?KQ$Mm ^OA Binding CAack^d LairUnate, PlcUd 

PdnCfRAH mCTA { INPUT ♦OUTPUT) 

peal MON (A) (4*At « ) *GiU*M) fDiU) *PT (A) 

PEAl. RI4) ,r(A) 

PEAL LP( 19) »OTAn<)) 

EQUIVALENCE (NON ,8) 

COMMON Kl»K2»K1tK4 
COMMON/AUX/H*P*PKl ♦PK2.0MO»X»y 
LPiU«o#n 
OTA n ) *0.0 
PEAO 2tK 1 fK2tK3,r>A 
2 E0PMAT(12) 

* K1 » Opoe« OF SYSTEM nF EQUATIONS 

* K? « MO. OF DISTINCT ‘fERNELS 

* K3 « NO. OF DATA POIMt$ 

* K4 « NO. OF DATA SETS TO BE EVALUATED 

* SET UP DATA POINTS 

AK*K3 

DO S Nsi,K3 

AN*N 

5 pT(N)=AN/AK 

* SET UP INTEORATION MArpiX 

M»K3-2 

N»K3-1 

AsKO 

Aal./(3.»A) 

no 10 K=?*M,2 
10 D(K)s2,*^a 

00 IS Ksl »N»2 
IS 0(K)=*+,«A 
0 (KD sA 

* calculate NONHOMOGENFOUS TERMS 

RHS=1 .0 
no 2? 1=1 »K2 
PRINT 9 
9 format (IH l) 

DO 999 lIsl.KA 
DO 35 N=l»K3 

35 NON(N) =RHS<»SQkT (PT (N) ) 

* CALCULATE KERNEL MATRICES 

call CONST (I) 

DO 20 N=1 »K3 
no ?0 M=1 .K3 

F(M.N.I)=FlJ(IfPTlM)»PT(N)) 

20 CONTINUE 

CALL CHANGE (FtG.D.T) 
call lined (G.B«C* K3) 

DO 40 L*1»K3 
PRINT 6fPT (L) »NnN(L) 

6 FORMAT (SX .F8.4.FIS.6) 

40 CONTI NUF 

LP 1 1 I>1 )=NOn(K3) 

0 T A ( I I ♦ I ) =P 
999 CONTINUE 

call LADINV (0TA,LH) 

22 CONTINUE 
END 


Mi'JCtlON n .A.m) 

COMMON/ AUX/H ♦ M ♦ PM « PK2 ♦ HMU » X t Y 
MXY?«2**lb 
nFL*0.2b«»(P-A» 
lF(nfL)A 0 »AS »50 
45 «;lMPaO.O 
PETURN 

50 CONTiNUf 

SA*zn«A)+ 7 l l»H) 

Sd*Z(I*At 2 .«» 0 FL) 

SCaZn ♦ A*l)FL) ♦? ( !♦ A>T,*nEL) 
SI*(DEL/3.)*(SA*2.<»SH»4.<»SC) 

IF(Sl.fO.O.O) GO Tn 45 
K»A 

35 SB»Sfl*SC 

DELs0.5*OEL 
^CxZ( 1 tA^OEL) 

J»K-1 

no 5 N=3»j.2 

AN*N 

5 SC*SC*Z< I*AtAN*f)EL) 

S2*(DEL/3.H^<5A42.»SB*4.»SC) 

niFsAssi is 2 -sn/si) 

FRsO.Ol 

TF(niF-FR) 30f25.25 
30 SlMPaS2 
PETIIPN 
25 K*?*K 

SUS2 

IF (K-MXYZ) ISt J5»40 

40 PRINT 42t 1 . A«H 

42 format (SXt** INT, DOES NOT CONVERGE «tl3»2F9.A) 
PRINT 60»X»Y 
ftO FORMAT(2FI0.5) 

DO 70 JxltlO 
OlPaJ 

DlPsDlP/10. 

VisZ(I»niP) 

PRINT 60«W 
70 CONTINUE 
CALL EXIT 
F.NO 


^'^1. PAOt II 
POOH QUAii,!; 


- 32 - 


SURROUTJNr UINtOU.BfTfN) 
R£At A(M»N) »P(N) »f <N) 
no 5 I»?.N 

s A( f .i)KAM*i)/An«n 
no 10 K»?«N 

M«K-1 

DO 15 t»l»M 

15 T(1)*A(I»K) 

no ?0 J*ltM 
A (JfK)»T( J) 

JlaJM 

no ?0 IsJlfN 

T ( 1 ) aJ ( I ) -A (I . J ) » A ( J ♦ K > 

20 CONTINUE 

A(K,K)sT<K) 

TF(K.E0.M go to 10 
M*K ♦ 1 

no 25 I*H»N 

25 A(I,K)aTm/A(K,K) 

10 CONTINUE 
• BACK SIJBSTlTUTe 
no 11 1=1 »N 
T n ) SB ( I) 

M* I ♦ I 

Ir (M.GT.M 1.-0 ro ?1 

no 10 jsM.N 

o( J)=B(J)-A(J*1)»T(1) 

10 CONTI Mur 
31 CONTINUE 
00 35 Isl ,N 
KsM+l-I 

R(K)=T (K)/A|KfK) 

Kl=K-l 

irCKl.EO.O) GO TO 15 
no 16 Jl=lfKl 
JsK-Jl 

T(J)=T(J>-A(J,K)i»HfK) 

36 rONTlNUt: 

35 CONTINUE 
oETURN 
END 
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FtJNCTIOM FU<I»A*H) 
CC-*HON/AI)X/H»P»PKl ♦PKa»RMU»X»Y 
X*A 
Y»H 

rF(4*9)5»lO»5 
10 FUaO.O 
PETIJRN 

5 SUMsSIMPd tO.OtS.O' 

FPsO.Ol 
DFL *5.0 
?0 UP*OEL*5.0 
* AOHLsSIMPi I *OEL.U?) 

OEL »UP 

TEST*AHS(ADDI./SIJM) 

SUMsSUM^AODL 
iFJTEST-fP) lS.20t?0 . 

1 ^ FUaSORT (X*Y) *SUM 
RETURN 
FNO 


SUPROUTINt CHANive (r,G.n. I ) 

peal F { 4 » 4, 1 ) ,(1 (A , 4 *) .0 ( 4 ) 

common K1 «KH,K3 .J' 4 

DO 10 Nal »K3 

DO 10 M=1 .K3 

r, (M,M) =r ( m.n, I ) <»o (N) 

10 CONTINUE 

DO PO Nal .Kl 
20 G(N.N)aG(N.N) ♦! .0 
RETURN 
FNO 
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rUMCTION BESJO(A) 
IFU-3.)5*5»10 
5 H=A»A/9. 

W*1 ,-^.2499997«H 
7*B*B 

WsVIti .P656?0«^Z 
7 = Z«P 

wsW-.3l63yft6»Z 

7*Z»B 

W-W* , 04U44 79<»Z 
2s7»f\ 

W=W-.003Q4u4»Z 

7*7«B 

PESJO=W*.0002l*»2 
RETIJPN 
10 fi=3./A 

Ws.797flP45A-.00000O77*B 
V = A-,78539«16-. 04 1^,6397*6 

WsW-.00b5274ttZ 
V=V-»000039S4*Z 
Z = 7»H 

w=w-,noon9‘=.i?«z 

V=V*.00?62573»2 
7 = Z»B 

Wsw».00l37?37«*Z 

VsV-.OOOSAiaS^Z 

Z-Z*»B 

W = W-.00072fl05<*Z 
V*V-.00029333»Z 
Z-Z»R 

WsW*.000l4476»7 

V*V^.00013S58«Z 

nESJOsW/SQPT (A) «COS(V) 

RETURN 

FNO 
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SUHPOUriNE CONSTdl 

COMMOli/ A«iX/H» P • PK 1 « PK * X ♦ Y 
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vC 

vr 




STATI«1TICS 

P&OGRAM LFf'f.^ j 

SfM LABFLPn COMMOM LFNbTH 7H 

47600P SCM USpO 


rUMrTtoN 2n«s) 

COMmON/AUX/H.P»KK 1 .PK2.HMU»X»y 
COMPLEX ACtALl fAL^.BAfSH 
complex G a « G0 ♦ B a ♦ HP t MC ♦ F ♦ 0 
Pls3.1A1592o 

ppsp«p 

AAsI . t7.»BMU*( 1 .-Pk1)/( I ,-PK2) 
nEOsAfl.* u ,-PKl ) / (Pp*h<*H*AA) 

SS = H»H<»AA/ {2,*»PI*PT<» « 1 ,-PKl )» ( 1 .♦BMli) ) 

XNllOs(PKI-PK2*AA«pw'2)/AA 

A0s(R4.SS)*OFO 

7Zs (P-SS ) * (P-SS) *0F 0»0E'0-A.»DE0 
G=CMPLX (Z2.0.0) 

ACaCSOPT(O) 

ALlsC.B^* (AB»AC) 

ALpsO.B* ( AR-AC) 

AL3s?.» (P*PEO*l n .-XNUO) 

SAsALE/ IP*DEO fl ,/SS) 

SBsALI/ (P^DEO^ l ./S«^) 

GArCSCjPT (G<>S»AL1 ) 

GBsCSOPT ( S^S*AL2) 

GCsSOPT (S‘»S + AL3) 

PAsP./n ,-XNU0<*XMUA)/(AL1-AL2) 
bB = GA»GA-XNUO<»S*S 
PC = C7B»GB-XI'JU(M»S*S 

FSHA*( ( 1 .-SA)«BH*RH/GA-(1 ,-SR)<»HC«HC/GP-2.»S*S*GC*(1 . 
II ?)/AL3) 

0=ocal(F) 

qA = aIMAG IF) - - 

jF(QA-0.0)«^f 10*5 

10 7= n-S) <»HESJO(S<‘X)»fiFS JO(S*Y) 

petdrn 

5 »PImT 9»P»S*F _ 

q FORMAT I 4F I 0.5) 

CALL EXIT 
FNO 
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-XNUO) (ALl-A 


i 


-rr-r!sr.t" 


^UBwOUT IMF LAPIMV (r.LAM^PHI ) 

C THIS PPOOHIAM FVAUMTKb Thfc' Cn£FK I C 1 M TS FOP SfPfFS 

r OF JACOHI OOLYNOMI&LS WHICH HtPHKStNTS A ( APLACE 

C INVFPSinN INTEGWAL 
PEAL HUL 

niMFMSION A(S0) tGLAM(bO) »PHl <S0) ♦C<A,S0) 

OIMFNSIDN HKUOl) .TT(IOI) 

COmmON/ 2/TI ,TF»nT fMM*9K»TT 
PE AO l»MN*MNtMM 

1 F0PMATni2) 

PEAO ?»TI*TF»OT 

2 FOPmAT OF10.5) 

PPIMT PP 

99 FOPMAT (1 HI) 

fall SPLICFIGLAMtPwI 
PHINT 10 I 

101 FOPMAT (/////5X«» GLAM PHI *) 

PPINT 1()2» (GLAMil) ♦PHI (I) ♦Ial»MM>, 

10? FOPMAT ISXtFlO.StSX.FlO.S) 

MllsMM-1 
PHIMT 99 

no 10 I s 1 ♦ MM 

PEAO 3,prT»DEL 

3 format (?F I O.S) 

PPIMT 9M*BFT»DEL 

98 FORMAT (V////5X»<»eETA s*F5.3»*> DELTA s«F5.3) 

no 11 l=i»mn 

AL=|. 

9=1 ./(AL^HET)/DFL 

CALL SPL INF (GLAM»PmI«Mm,C*S»G) 

F=r,»s 

lFUL-?.)Hl»fl?«63 
81 A(1)=(1.*BET)«0EL*F 

go to 'll 

e? A (?) = ( (2.>HET)*r)EL<^F-A(l) )»(3.>HET) 

GO TO 11 
• 83 continue 
tops], 

L1=L-1 
ALl sLl 

no I? J=1 fLl 
AJ=J 

top=aj*top 

1? continue 
L2=?*L-1 
hOTsI. 

no 13 J=L»L2 
A J = J 

HUT=(AJtRET) *fiOT 
13 rOMTINUF 
mul=bot/tod 
SUM=0.0 
no 14 N=1 •Ll 

ANsrj 

IF ( AN-?. ) 85 f 86 » 87 
85 TOn=l. 




AA TOO*AU 
GO TO AA 
A7 CONTI NUF 
TOOsl. 

no 15 J=ICM,L1 

A J3 J 

TOOsAJoTOO 
15 COMTINUF 
bA CONTiNlIh 
POOsl. 

JAsLl*N 
no 16 J=L#JA 
AJsJ 

BOOsBOO*(AJ^BET) 

, 16 CONTINUE 
COsTOO/BOO 
SUM=SUM^CO<»A(N) 

14 CONTINUE 

A(L)®MUL*(nEL*F-SUM) 

11 CONTINUE 

CALL JACSER(OEL»A«f-ET) 

10 CONTINUE 
49Q CONTINUE 
RETURN 
END 

5UHB0UTINE JACSEP (n»C»H)‘" 
dimension C(50) »Sr(50) »P (50) 
dimension PK(lOl) fTT(lOl) 

COMMON/2/TI*TFtDT*MN»BK.TT 
TT ( I )s0.0 
PK(J )s0.0 
LMsl 
T = TI 

12 T=T*OT 
X*2.»EXP(-0«T)-1. 

CALL JACOBI (MN*X*BtP) 

5Fn)sC(l )*P(1) 
no 10 l= 2 *mn 
LUL-l 
AL=L 

• SF(L)=SF(LI)>C(L)*oU.) 

10 continue 

LM=lMM 
PK(LM)=5F(5) 

TT (LM)=T 

IF(T.LE.TF) go to 12 
print 97 

97 FORMAT (///// 5X »<» T K T K T K 

IT K *) • . 

DO 31 MY=1 »25 
MA=MY»1 

M0=MA^?5 

mCsmP*2S 

mO=MC>25 

PRINT 9B, TT (MA) .BK (MA) f TT (MB) »BK (MH) »TT (MC) »BK (MC) »TT (MD) fPiK (MD) 

96 FORMAT (5X f F5.2t1X»r7.5» 3X»F5.2»3X»F7.5»3X»F5.?»3X»F7.5 v3X*F5.2,3X. 

1F7.S) 

31 CONTINUE 
RETURN 
FNO 


- 39 - 


o o 


SUHPOOTINiE JACOWI 

THIS PHOOHAH CALCKI ATES JACOH I POLYNOMIALS OF OROEW 
K-1 WITH APG X AND PAHAMETER H GT -1 
dimension Pb(N) 

ANsN 

IF<an- 2.) I 
I P«(l)=l. 

RETURN 
? PBUIsl . 

PB( 2 )=X*e<*(l.-X»/ 2 , 

RETURN 

3 PS0sR*8 
RONFsB*!., 

PBU) = l. 

PB(2)=X-R*(l.-X)/2. 

DO A K=3fN 
AK=K 

AKLsAK-l. 

AK?=AK-2. 

Kl=K-l 

K2=K-2 

COU< (2.*AKl) tB)<*X ■ 

COU( {3.«AK2) ♦B)*COl 

COls( i2,*AK2) ♦BONE) *(C01-8SQ) 

C02a2.<»AK2*(AK2>e)»( (2.*AKl)^b) 

COs?.*AKl<» (AKl*R)*( (2.*AK2) *8) 

4 Pe(KJ = (C0l*Pfl(Kl)-C02*PH(K2M/C0 
RETURN 

FNO 


Oil 
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SUHROUT 1 ME SPL I ME (X f Y • M ♦ C ♦ A I M f ♦ Y 1 NT ) 

DIMENSION X(S0) *Y(SO> .C(4»50) 

1F(XINT-Xm)l00»ll 

10 YiNTsYd) 

RETURN 

1 1 CONTINUE 
TF(X<M)-XIMT)1»12«|3 

12 YlMTsY(M) 

RETURN 

13 CONTINUE 
KaM/2 

N = M 

? continue 

IFU(K)-X1NT)3*U»S 

14 yInT=Y{K) 

RETURN 

3 CONTINUE 

IF(XINT«X{KM) )4»1«»7 

15 YINT=Y{K^1) 

RETURN 

4 CONTINUE 

YlMTs(X IK^l )-XlNT)'MC( 1 . K ) » ( X ( K ♦ 1 ) -X INT ) ••2^G H ♦ K ) ) 
YiNTaYINT* (X1NT-X(<) ) * I C ( 2 ♦ K ) I X I NT-X iK ) )<»*2*C(4«K) ) 
PETllPN 

s continue 

IF(X (K-n -XINT lh*17 

6 KsK-l 
GO TO 4 

16 YlNTaY(K-l) 

RETURN 

17 N=K 
KaK/2 
GO TO 2 

7 I.LsK 
K=(N+K)/2 

8 CONTINUE 

lF(X<K)-XINT)3»14n8 

18 continue 

IF(X(K-1)-XINT)6»16»19 
t9 N*K 

Ka(Ll.*K)/2 
GO TC 8 

1 print 101 

101 format (« OUT OF RANGE FOR INTERPOLATION ») 

STOP 

END 
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SUPROUTIME SPLlCt:(^«Y*M»C) 

DIf-'FNSlON X(50) fYC^O) tO^SO) fPiSO) *e (50) *CU*50) 
niMENSION A(50»3) ♦P(50) fZiSO) 

no ? Ksl*MM 
0(K)sX(K*l )-X(K) 

P(t<)rO(K)/6. 

? r(K)s(Y(K^l)-Y(K) )/0(K) 
no 3 Ka2«MM 
3 «(K)aE(K)-E(K-l) 

An,2)=-i.-o(i)/i:)(?) 

A(i*3)=rHn/o(2) 

A(p,3)sp(2)-pn)^An»3) 

A(?.2)=2.‘»(P(1)*P(P) )-P(l)‘»A(U2) 
A(?.3)=A<2»3)/A(2,?) 

R(2)=H(?)/A(2.2) 

DO A Ks3»MM 

A(K,2)=2.«(P(K-1)+P(K) )-P(K-l)*A(K-lf3) 
R<K)=B(K)-P(K-1 )*8(K-l) 

A(K,3)=P(K)/A(K.2) 

A y(K)=B(K)/A(Kf2) 

OsO(M-2)/()(M-l ) 

A (M, I ) =1 . ♦n«'A (M-20) 

A <m, 2) s-O-A (M, 1 ) «A (M-l «3) 

B (M> =B (M-2) -A (M, 1 )<^e (M-l ) 

7(M)sB(M)/A(Mf2) 

MN=M-2 
no 6 IsltMN 

KaM-I 

(S 7(K)sH(k)-a(K,3)<»Z(K*1) 

7(1 )a-A( 1 f?)»Z(2)-A(l t3)*Z(3) 

DO 7 K=1 tMM 
0*1 ./(6.*0(K) ) 

C(l»*<)=Z(K)^0 

C(?»K)=7 (K*l ) <»0 

C (3»K)=Y (K)/D(K)-7(K) *P(K) 

7 C(4»K)sY(K*1)/D(K)-Z(K>1 )*P(K) 

PETURN 

END 


IntziUitif Factor 


MUiJ/MUl * 50»00 NiM ■ ,30 MU2 * ,30 
A/H «1.00 C21/PA ■ .02 


• 2500 .0338f=.5 

.5000 ,0598?9 

.7500 .090175 

1.0000 .257758 


MU2/MU1 = 50.00 NUl = ,30 NU2 * ,30 
A/H =1,00 C21/PA = ,04 


.2500 

.5000 

.7500 

1.0000 


.069541 
. 119396 
.10301 3 
.352715 


MU2/MU1 = 50.00 NUl - .30 NU2 = .30 
A/H =1,00 C21/PA = .06 


.2500 .103272 

•5000 .169658 

•7500 .249217 

1.0000 .412993 
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t/e/Uvcition Equation (25). Equation (25) can be derived by first expressing 
equation (24) In the form 


p„h3 52 j^_ 




* I « // It «<^Jy + To * 

'o 


(59) 


In which dW/3t can be written as 


+ + _5L!l!!2 


dr 


ar. 


at dr^ at afy at 


(Mv It + H,, + Qv) 


ary, at 




X ax xy ay ^x^ at 


^^xy ax ^y ay ^y^ at^ 


^<^xlxMy|7)|| 


(60) 


Denoting n and s as the normal and tangential direction, equation (60) may be 
integrated to yield 


// f dxdy = ^ . II Q,)ds - // (|;i . ■ Q,) 


di(i 


di/i 


BK 3H, 


+ ( ^*^xy . n ) + HifHv 

* JT ^Jy’ V ^W‘ * ay 


(6U 


Putting equation (61) into (59) and observing the relations in equations (15), 
the expression for T+? in equation (25) is obtained. 
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